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1.  INTRODUCTION 


During  the  past  -fifteen  years,  a  great  deal  o-f  research  has  -focused 
upon  the  task  of  devising  closure  approximations  for  the  long-time- 
aver aged  Navi bt — Stokes  equations  suitable  for  predicting  properties 
of  turbulent  flows.  Prior  to  1968,  virtually  all  turbulence  clos¬ 
ure  schemes  were  " incomplete",  i.e. ,  their  implementation  required 
some  advance  knowledge  about  the  flowfield  under  consideration  in 
order  to  obtain  a  solution.  The  best  known  incomplete  turbulence 
model  is  the  mixing-length  model  (Ref  1).  This  model  is  incomplete 
as  the  appropriate  form  of  the  mixing  length  must  be  determined  em¬ 
pirically  for  each  new  application;  in  general,  it  cannot  be  speci¬ 
fied  a  priori. 

In  1968,  the  first  "Stanford  Olympiad"  (Ref  2)  was  held  to  test  ex¬ 
isting  turbulence  models  against  the  best  experimental  data  avail¬ 
able.  The  data  base  was  confined  to  incompressible  two-dimensional 
boundary  layers.  Interestingly,  the  competition  was  more-or-less 
won  by  the  "complete"  model  of  Bradshaw,  et  al  (Ref  3). 

Perhaps  spurred  on  by  the  success  of  the  Bradshaw  model ,  the  trend 
in  turbulence  modeling  since  the  first  Stanford  Olympiad  has  been 
toward  development  of  complete  models.  For  clarity,  note  that  the 
terminology  "complete  model  of  turbulence"  as  used  in  this  report 
means  a  set  of  equations  which  can  be  used  to  predict  a  given  turb¬ 
ulent  flow  with  no  advance  information  other  than  boundary  condi¬ 
tions  required  in  order  to  achieve  a  solution.  The  terminology  is 
not  intended  to  imply  anything  with  regard  to  the  range  of  applica¬ 
bility  of  the  theory. 

Over  the  past  decade,  the  most  vigorous  modeling  efforts  have  been 
conducted  by  Donaldson,  et  al  (Ref  4),  Launder,  et  al  (Ref  5),  and 
Wilcox,  et  al  (Refs  6-9).  Recognizing  the  substantial  progress  the 
various  researchers  seemed  to  be  making,  the  second  "Stanford  Olym¬ 
piad"  was  held  in  1980  and  1981  (Ref  10).  This  time,  however,  the 
scope  of  the  experimental  data  was  expanded  tremendously  to  include 
complicating  effects  of  compressibility,  streamline  curvature,  sur — 
face  mass  transfer,  boundar y-1 ayer  separation,  secondary  motions, 
etc.  That  is,  virtually  every  complicating  effect  known  to  man  was 
included  in  the  Olympiad  if  experimental  data  of  reliable  quality 
existed. 

From  this  researcher’s  viewpoint,  the  results  of  the  Second  of  the 
two  Stanford  Olympiads  were  at  once  very  encouraging  and  dissapoin- 
ting.  On  the  one  hand,  the  state  of  the  art  has  been  shown  to  have 
advanced  far  beyond  the  wildest  dreams  of  evaluators  of  the  First 
Olympiad.  It  was  hard  to  imagine  in  1968  that  separated  flowfields 
could  be  routinely  predicted  with  any  degree  of  accuracy  Just  fif¬ 
teen  short  years  later.  COf  course,  turbulence  modelers  should  re¬ 
ceive  only  part  of  the  credit;  magnificent  advances  in  numerical 


methods  such  as  those  of  MacCormack  (Ref  11)  have  played  a  very  im¬ 
portant  role  to  say  the  least!]  On  the  other  hand,  while  such  pre¬ 
dictions  can  be  routinely  made,  obtaining  results  consistent  with 
experimental  measurements  is  not  nearly  as  routine.  Far  worse,  it 
is  not  even  clear  from  the  results  presented  at  Stanford  Olympiad  2 
that  we  can  predict  effects  of  an  adverse  pressure  gradient  on  the 
turbulent  boundary  layer  any  more  accurately  than  we  could  fifteen 
years  ago.  Clearly,  our  progress  in  turbulence  modeling  has  been 
a  bit  uneven. 

In  light  of  this  situation,  this  study  has  been  initiated  with  the 
objective  of  first  taking  a  modest  step  backward  by  reviewing  and 
assessing  the  original  closure  approximations  for  the  class  of  tur — 
bulence  models  known  as  two-equation  models,  viz,  closure  being  ac¬ 
complished  using  the  long-time-averaged  Navier — Stokes  equations  and 
two  additional  differential  equations.  The  rationale  for  starting 
at  what  would  seem  to  be  a  very  elementary  level  stems  from  a  key 
observation  made  at  the  Second  Olympiad,  viz,  the  greatest  amount 
of  uncertainty  and  controversy  over  two-equation  and  higher — order 
models  lies  in  the  scale-determining  equation.  It  is  even  unclear 
as  to  what  the  optimum  choice  of  dependent  variables  is  for  a  two- 
equation  model.  As  a  result  of  this  study,  we  feel  we  have  found 
the  optimum  choice  and,  based  upon  this  choice,  we  have  postulated 
a  new  two-equation  turbulence  model. 

Section  2  summarizes  the  new  model,  including  arguments  which  set 
values  of  all  but  two  of  the  closure  coefficients  appearing  in  the 
postulated  equations.  Section  3  presents  results  of  a  perturbation 
analysis  of  the  incompressible  defect  layer,  including  effects  of 
pressure  gradient.  Predictions  of  the  new  model  are  compared  with 
those  of  the  Jones-Launder  (Ref  5)  and  the  Wilcox-Rubesin  (Ref  9) 
models.  Section  4  uses  perturbation  methods  to  analyze  the  viscous 
sublayer,  including  effects  of  surface  roughness  and  surface  mass 
injection.  Section  5  presents  results  of  application  of  the  model 
to  five  incompressible  boundary  layers.  The  concluding  section  in¬ 
cludes  a  summary  of  results  obtained  in  this  study. 


2.  EQUATIONS  OF  MOTION 


In  this  section,  w  stats  ths  postulated  squat  ions  o-f  sot  ion,  in¬ 
cluding  sstablishsd  valuss  of  all  closure  coefficients.  Physical 
interpretations  of  turbulence  field  properties  are  given  and,  addi¬ 
tionally,  arguments  are  presented  Mhich  have  been  used  in  setting 
values  of  several  of  the  closure  coefficients. 


2.1  POSTULATED  EQUATIONS 

For  general  incompressible  turbulent  fluid  flows,  the  turbulence  mo¬ 
del  equations  are  as  follows. 


Mass  Conservation 


Momentum  Conservation 


+  u, 


Turbulent  Mixing  Energy 


1*  +  u  It 
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Turbulent  Dissipation  Rate 
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where  t  is  time  and  x^  is  the  position  vector*  u ^  denotes  long-tim 
averaged  velocity  vector*  p,  p  ,  and  v  are  mean  pressure,  density 
and  kinematic  viscosity*  and  is  the  Reynol ds-stress  tensor.  Th 
turbulent  mixing  energy,  k,  and  the  turbulent  dissipation  rate,  u  , 
are  needed  to  define  the  eddy  diffusivity,  v,p  ,  which  is  given  by 

V-,  ■  y*  k  /  w  CS) 


m  strain 


The  Reynol ds— stress  tensor  is  assueed  proportional  to  the  eea 
rate  tensor,  S^j  ,  so  that  we  write 


(6) 


where,  by  definition. 


if  iii 


+ 
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Several  closure  coefficients,  viz,  3,  3*,  Y,  y*,  a,  and  a*  appear 
in  Equations  <3— 3) .  A  key  objective  of  this  study  has  been  to  review 
typical  arguments  used  in  establishing  values  of  such  coefficients  in 
a  eodel  of  this  type.  In  the  next  subsection  and  in  later  sections, 
the  argueents  are  presented.  For  convenience,  the  values  are  summar- 
ized  in  the  following  equations. 


6  =  3/40  ,  3*  =  9/100 

Y  »  3/9  ,  Y*  -  1 

a  *  1/2  ,  a*  »  1/2 


Before  proceeding  to  further  discussion  of  the  closure  coefficients, 
it  is  worthwhile  to  pause  and  discuss  the  fora  of  the  model  equations 
and  the  physical  meanings  of  the  quantities  k  and  u>.  As  in  other 
two-equation  models  of  turbulence,  the  quantity  k  represents  a  mea¬ 
sure  of  the  kinetic  energy  of  the  turbulence.  Whether  k  is  specifi¬ 
cally  identified  as  being  the  exact  kinetic  energy  of  the  turbulence 
or  alternati vely  as  the  kinetic  energy  of  the  fluctuations  in  the  di¬ 
rection  of  shear  (Ref  8)  is  not  critically  important.  All  we  require 
on  physical  grounds  is  that  k  be  proportional  to  the  square  of  the 
velocity  at  which  local  turbulent  mixing  occurs.  The  second  quantity 
introduced  in  the  model,  <*>,  is  referred  to  as  the  turbulent  dissipa¬ 
tion  rate.  Its  dimensions  are  inversely  proportional  to  time  and  it 
is,  in  fact,  the  same  variable  used  by  this  author  in  all  prior  turb¬ 
ulence  modeling  studies.  Perhaps  the  simplest  physical  interpretation 
of  u  is  that  it  is  the  ratio  of  the  turbulent  dissipation,  e,  to  the 
turbulent  mixing  energy.  A1 ternatlvely,  «  is  the  rate  of  dissipation 
of  turbulence  per  unit  energy. 

As  is  obvious  from  inspection  of  Equation  (3),  the  equation  for  k  is 
modeled  directly  after  the  exact,  long-time-averaged  equation  for  the 
turbulent  kinetic  energy.  On  this  point,  the  model  is  consistent  with 
virtually  all  other  two-equation  models.  The  second  of  the  two  model 
equations  is  similar  in  form  to  the  equation  for  k.  Although  it  adds 
no  rigor  to  the  approach,  the  equation  for  w  can  be  regarded  as  the 


■odtltd  form  of  thm  mquation  which  would  result  from  (a)  writing  thm 
exact  aquations  for  turbulant  kinetic  anmrgy  and  dissipation  and,  <b) 
making  thm  formal  change  of  dependent  variables  defined  by 


ui  *  e/(B*k) 


(9) 


The  primary  difference  between  the  model  postulated  in  this  study  and 
those  of  this  author's  prior  research  is  the  form  of  the  equation  for 
w.  Host  notably,  past  studies  have  written  the  equation  in  terms 
of  the  square  of  u.  Interestingly,  the  first  two-equation  model  in 
which  the  variables  k  and  u  were  used  was  postulated  by  Kolmogorov 
<^ef  12)  and  his  equation  for  u  was  written  in  terms  of  u  rather  than 
u  •  The  reason  for  this  choice  will  become  quite  clear  in  Section 
3  where  we  analyze  model -predicted  structure  of  the  defect  layer. 


2.2  ESTABLISHING  CLOSURE  COEFFICIENT  VALUES 


In  this  subsection,  we  present  straightforward  arguments  from  which 
values  of  the  four  closure  coefficients,  8  ,  8*,  y  ,  and  y*  can  be 
established.  A  review  of  the  arguments  generally  presented  indicates 
the  following  are  as  physically  sound  as  possible  within  the  context 
of  two— equation  turbulence  models. 

Considering  first  the  coefficient  y*f  we  can  rewrite  Equations  (3-6) 
in  terms  of  the  quantity  u/y*  .  Inspection  of  the  resulting  equa¬ 
tions  shows  that  this  rescaling  of  u  is  equivalent  to  setting  y*»  1. 
Hence,  with  no  loss  of  generality,  we  conclude  that  the  value  of  y*  is 
indeed  unity. 

Next,  we  turn  to  the  ratio  of  8*  to  8  •  For  decaying  homogeneous , 
isotropic  turbulence.  Equations  (3-4)  simplify  to 


dk/dt  -  -  6*ku> 

du/dt  ■  -  8u2 


(10) 


from  which  the  asymptotic  solution  for  k  is  readily  found  to  be 


k  t”6#/B  (11) 

-0 

Experimental  observations  indicate  that  k  <v  t  for  decaying  homo¬ 
geneous,  isotropic  turbulence  which  implies  0*/g  “  b/3  which  is  con¬ 
sistent  with  Equations  (8). 

Values  for  the  coefficients  y  and  8*  can  be  established  by  examining 
the  so-called  "wall  layer”.  The  wall  layer  is  defined  as  the  portion 
of  the  br  -ndary  *  jyer  sufficiently  distant  from  the  surface  that  mol¬ 
ecular  vl  ■  ii  is  negligible  relative  to  eddy  viscosity,  yet  close 
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enough  for  convective  effects  to  be  negligible  relative  to  the  rate 
at  which  the  turbulence  is  being  created  and  destroyed.  In  the  1.^ sit¬ 
ing  case  of  a  constant-pressure  boundary  layer.  Equations  (2-7)  sim¬ 
plify  to 

0  *  fy  [VT  fy] 

0  *  8**“  +  a*fy  [VT  I?] 

0  .  y[|^]  -  S  +  °  |y[vT  |f] 

Me  seek  the  conditions  under  which  these  simplified  equations  yield 
a  solution  consistent  with  the  law  of  the  wall,  i.e. ,  velocity  vary¬ 
ing  linearly  with  the  logarithm  of  distance  from  the  surface.  As  can 
be  easily  verified.  Equations  (12)  do  indeed  possess  a  solution  con¬ 
sistent  with  the  law  of  the  wall,  viz, 

UT 

u  =  —  Hn(uTy/v)  +  constant  \ 

k  =  u*  //&*  |  03) 

a)  =  uT  /(/£*  k  y)  ' 

where  uT  is  the  conventional  frictional  velocity  and  k  is  Karman’s 
constant.  There  is  one  constraint  imposed  in  the  solution  to  Equa¬ 
tions  (12),  namely,  a  unique  relation  exists  between  the  implied  val¬ 
ue  of  Karman’s  constant  and  the  various  closure  coefficients.  Speci¬ 
fically,  the  following  equation  must  hold. 


y  *  3/e*  -  aicV/e*  (i4) 


Additionally,  note  that  the  Reynolds  shear  stress,  r  »  is  constant  in 
the  wall  layer  and  is  equal  to  uf  .  Finally,  inspection  of  Equations 
(13)  shows  that  this  implies  T/ka/g*  in  the  wall  layer.  A  variety 
of  experimental  measurements  (Ref  13)  indicate  the  ratio  of  t  to  k  is 
about  3/10  in  the  wall  layer.  Thus,  the  predicted  wall-layer  solution 
is  consistent  with  experimental  observations  provided  3*  «  9/100. 

In  summary,  the  arguments  presented  in  this  subsection  are  sufficient 
to  uniquely  set  the  values  of  3  »  3**  end  y*.  Also,  Equation  (14) 
determines  Y  in  terms  of  the,  as  yet  undetermined,  value  of  a  .  As 
a  byproduct  of  analysis  in  the  next  section,  values  for  the  two  clo¬ 
sure  coefficients  a  and  a*  will  be  established. 


(12) 


3.  DEFECT-LAYER  ANALYSIS 


In  this  section  wb  use  singular  perturbation  methods  to  analyze  mo¬ 
del-predicted  structure  of  the  classical  defect  layer.  The  analysis 
presented  is  a  general ization  of  that  done  by  Wilcox  and  Traci  (Ref 
8).  In  contrast  to  the  Wilcox  and  Traci  analysis,  effects  of  pres¬ 
sure  gradient  have  been  included.  Additionally,  the  analysis  has 
been  done  for  three  turbulence  models,  vizi  the  model  postulated  in 
Equations  (1-8);  the  Wi 1 cox -Rubes in  (Ref  9)  model;  and  the  Jones— 
Launder  (Ref  5)  model.  First,  we  review  details  of  the  perturbation 
solution  procedure.  Next,  we  compare  solutions  for  the  three  models 
in  the  absence  of  pressure  gradient.  Then,  effects  of  pressure  gra¬ 
dient  are  studied  for  the  three  models.  Finally,  we  Justify  the 
values  chosen  for  the  closure  coefficients  o  and  cr*  . 


3.1  PERTURBATION  SOLUTION 

In  the  past,  the  only  detailed  analyses  of  the  defect  layer  for  any 
turbulence  model  have  been  those  of  Bush  and  Fendell  (Ref  14  -  for 
the  mixing-length  model)  and  Wilcox  and  Traci.  In  neither  case  were 
effects  of  pressure  gradient  delineated.  In  this  section  we  extend 
the  Wilcox-Traci  analysis  to  include  pressure  gradient. 

To  study  the  defect  layer,  we  seek  a  perturbation  solution.  The  ex¬ 
pansion  proceeds  in  terms  of  the  ratio  of  friction  velocity  to  the 
boundary-layer — edge  velocity,  uT/Ue,  and  the  dimensionless  vertical 
coordinate,  r|  ,  defined  by 

u  y 

n  *  TF  <15> 

e 

where  6*  is  displacement  thickness.  For  the  sake  of  brevity,  we 
confine  details  of  the  expansion  procedure  to  Appendix  A.  It  is  in¬ 
structive  to  note  that  the  velocity  is  given  by 


(n)  + 


(16) 


which,  to  order  u  /U 

t  e 


can  be  rewritten  as 


U 


e 

u 


-u 


T 


f(y/&)  ;  A  ■  U  fi*/u 

e  x 


( 17> 


The  coordinates  appearing  in  Equation  (17)  are  the  classical  de¬ 
fect-layer  coordinates.  Additionally,  it  is  important  to  note  that 


pressure  gradient  appaars  in  the  aquations  of  motion  in  dimansion- 
lass  -form  as 


&T  *  6*(dp/dx> /  tw 


(18) 


where  tw  is  the  surface  shear  stress.  Colas  and  Hirst  (Ref  2)  refer 
to  3^,  as  the  equilibrium  parameter. 

In  order  to  solve  the  defect-layer  equations,  we  have  used  an  im¬ 
proved  version  of  the  implicit  time  marching  program  developed  by 
Wilcox  and  Traci.  That  is,  we  add  unsteady  terms  to  each  of  the 
equations  of  motion,  make  an  educated  guess  of  the  solution  and  in¬ 
tegrate  over  time  until  the  solution  displays  negligible  temporal 
variation.  Computations  have  been  done  using  31,  SI  and  71  points 
in  the  grid  normal  to  the  surface  to  determine  solution  sensitivity 
to  the  numerical  algorithm.  For  all  three  models  considered,  dif¬ 
ferences  between  solutions  obtained  using  31  and  71. mesh  points  are 
less  than  two  percent.  Results  of  all  computations  given  in  this 
section  have  been  obtained  using  51  mesh  points. 


3.2  FLAT-PLATE  BOUNDARY  LAYER 

Figure  1  compares  numerical  predictions  of  the  three  models  with 
corresponding  experimental  data  of  Weighardt  (Ref  2).  (Note  that  in 
the  new  model  computation  we  use  a  »  o*  =  1/2$  we  defer  any  fur — 
ther  discussion  of  the  appropriate  values  to  Subsection  3.4. >  The 
experimental  data  presented  are  those  at  the  highest  Reynolds  num¬ 
ber  for  which  data  are  reported.  This  is  consistent  with  the  de¬ 
fect-layer  solution  which  is  formally  valid  for  very  large  Reynolds 
number.  Numerical  results  are  shown  for  three  models,  viz,  the  new 
(NEW)  model,  the  Wilcox-Rubesin  (W-R)  model,  and  the  Jones -Launder 
( J-L)  model . 

As  shown,  all  three  models  predict  velocity  profiles  which  differ 
from  measured  values  by  no  more  than  about  three  percent  of  scale. 
Interestingly,  the  new  model  shows  the  smallest  differences  from 
the  Weighardt  data.  Additionally,  skin  friction,  c^,  can  be  infer¬ 
red  from  the  defect-layer  solution  (see  Appendix  A).  Corresponding 
computed  and  measured  values  are  summarized  in  the  insert  on  Figure 
1|  the  largest  difference  is  less  than  three  percent.  Thus,  based 
on  analysis  of  the  constant-pressure  defect  layer,  there  is  little 
difference  amongst  the  three  models. 


3.3  EFFECTS  OF  PRESSURE  GRADIENT 

Turning  now  to  the  effect  of  pressure  gradient,  we  have  computed 
defect-layer  solutions  for  the  equilibrium  parameter,  gip  ,  ranging 
from  -0.5  to  >9.0.  Note  that  positive  $ij  corresponds  to  an  adverse 
pressure  gradient.  The  choice  of  this  range  of  Srp  has  been  dictated 


(U  -u)/u 


Figure  1.  Comparison  of  computed  and  measured  velocity 

profiles  in  defect-layer  coordinates  for  a  con¬ 
stant-pressure  boundary  layer. 


by  the  rcquirnent  of  the  perturbation  solution  thatgip  be  constant. 
This  is  as  wide  a  range  as  we  have  been  able  to  find  for  Mhich  ex- 
peri  eental  data  have  been  taken  with  6^  eore-or-less  constant. 


Figure  2  coepares  computed  wake  strength,  if  ,  with  values  inferred 
by  Colts  and  Hirst  from  experimental  data.  For  the  sake  of  clarity, 
note  that  the  wake  strength  appears  in  Cole’s  composite  law-of-the- 
wall/wake  profile,  viz. 


B  ♦  ■  ¥  sin2[r  f] 


(19) 


Appendix  A  summarizes  the  method  used  to  infer  if  from  the  numeri¬ 
cal  predictions. 

Inspection  of  Figure  2  reveals  provocative  differences  amongst  the 
three  models.  Most  notably,  the  new  model  yields  wake  strengths 
closest  to  values  inferred  from  data  over  the  complete  range  consi¬ 
dered.  Consistent  with  predictions  of  Chambers  and  Wilcox  (Ref  IS), 
the  Jones -Launder  model  exhibits  the  largest  differences,  with  pre¬ 
dicted  wake  strength  502-100%  lower  than  inferred  values  when  BT  is 
as  small  as  two! 

Figure  3  compares  computed  velocity  profiles  for  *  9  with  exper¬ 
imental  data  of  Clauser  (Ref  2)  (or  &T  *  B.7.  As  with  the  constant 
pressure  case,  computed  and  measured  skin  friction  are  included  in 
the  insert.  Consistent  with  the  wake-strength  predictions,  the  new 
model  yields  a  velocity  profile  and  skin  friction  closest  to  meas¬ 
urements  while  the  Jones-Launder  model  shows  the  greatest  differ¬ 
ences.  The  Wilcox-Rubesin  profile  and  skin  friction  lie  about  mid¬ 
way  between  those  of  the  other  two  eodels.  Using  Clauser’s  value 
for  displacement-thickness  Reynolds  number,  the  velocity  profiles 
are  replottod  in  sublayer  coordinates  in  Figure  4.  As  shown,  dif¬ 
ferences  amongst  the  models  are  as  pronounced  in  sublayer  coordi¬ 
nates  as  they  are  in  defect-layer  coordinates. 

The  explanation  of  the  Jones-Launder  model’s  poor  performance  for 
adverse  pressure  gradient  can  be  developed  from  inspection  of  the 
asymptotic  behavior  of  solutions  as  n  ->  0.  For  the  three  models 
tested,  the  velocity  behaves  as  follows. 

Vu  i 

— —  -  —  Ann  +  A  -  8^,Bn  Ann  +  ...  120) 


where  the  constants  A  and  B  are  summarized  in  Table  1.  Note  that, 
while  the  coefficient  A  is  determined  as  part  of  the  solution  (from 
the  integral  constraint  that  mass  be  conserved),  the  coefficient  B 
follows  directly  from  the  limiting  form  of  the  solution  as  n  °* 
As  seen  from  Table  1,  B  is  largest  for  the  Jones-Launder  model  and 
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Figure  2.  Comparison  of  computed  and  measured  wake  strength  for  equilibrium 
boundary  layers  with  pressure  gradient. 


MODEL  I 103Cf 


Table  1. 


Summary  of  Coeff i dents  A, 
in  Equations  <20-21 )  -for 


B  and  L 
St-9. 


Model 

A 

B 

L 

New 

13.1 

10.6 

-19.8 

Wi 1 c ox -Rubesi n 

9.8 

23.4 

-32.6 

Jones-Launder 

5.4 

54.8 

-61.1 

is  smallest  -for  the  New  model.  The  presence  o-f  the  r)  Anri  term 
gives  rise  to  the  in-flection  in  the  velocity  profile  as  rj  ->  0  that 
is  most  pronounced  for  the  Jones -Launder  model.  In  terms  of  turbu¬ 
lence  properties,  the  turbulence  length  scale,  l  *  kV2/u>  = 
behaves  according  to: 


l 


-> 


+ 


8TLn£nn  + 


(21) 


where  the  coefficient  L  is  summarized  in  Table  1.  Again,  we  see 
that  the  contribution  of  the  n  Inn  term  is  largest  for  the  Jones- 
Launder  model  and  smallest  for  the  New  model.  Thus,  in  the  presence 
of  adverse  pressure  gradient,  the  JonesHLaunder  1  tends  to  be  too 
large  in  the  near-wall  region.  Note,  of  course,  that  this  short¬ 
coming  is  not  evident  in  the  constant-pressure  case  which  has  8^  O. 

The  manner  in  which  the  New  model  achieves  smaller  values  of  £  than 
does  the  Jones-Launder  model  can  be  seen  by  changing  dependent  var¬ 
iables.  That  is,  starting  from  the  k  -  w  formulation  and  defining 
e  *  8*k(i)  ,  we  can  deduce  the  following  equation  for  e  implied  by 
the  New  model . 


de 

dt 


<  1+Y>  k 


<l+8/B»>  £2 


2cv  ife 
*  vT3y  3y 


(22) 


All  but  the  last  term  on  the  right-hand  side  of  Equation  (22)  are 
identical  in  fora  to  those  of  the  Jones-Launder  model  (see  Appendix 
A).  The  last  term  is  negligibly  small  as  n  ->  0  for  constant-pres¬ 
sure  boundary  layers  because  k  ->  constant  as  n  *  y/A  ->  0.  How¬ 
ever,  3k/3y  is  nonvanishing  when  0  and  9(e/k)/3y  generally  is 
quite  large  as  n  “>  0.  The  net  effect  of  this  additional  tern  is 
to  suppress  the  rate  of  increase  of  l  close  to  the  surface. 

As  a  final  comment,  note  that  all  of  our  computations  have  usmd  the 
model -predicted  behavior  (e.g..  Equations  20-21)  as  "wall -function” 
type  boundary  conditions  for  rt  ->  0.  Using  other  empirical  wall 
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■functions  presumably  would  improve  Jones— Launder  predictions.  How¬ 
ever ,  the  asymptotic  behavior  (e.g.,  inflected  velocity  profile) 
inherent  to  the  model  ultimately  must  prevail  at  high  Reynolds  num¬ 
ber  if  the  point +of  application  of  the  wall  functions  remains  con¬ 
stant  atf  say,  y  '■  80.  To  understand  this  point,  one  need  only  note 
that,  by  definition,  n  is  related  to  y+  by* 

H  -  y+  /Re6.  (23) 


Hence,  suppressing  the  asymptotic  behavior  inherent  to  the  model 
would  require  using  wall  functions  to  increasingly  larger  values  of 
y+  as  Reynolds  number  increases. 


3.4  ESTABLISHING  CLOSURE  COEFFICIENT  VALUES 

Unlike  the  four  closure  coefficients  discussed  in  Subsection  2.2, 
we  have  been  unable  to  find  satisfactory  arguments  to  establish  the 
values  of  a  and  a*  prior  to  performing  any  numerical  computations. 
However,  we  have  found  from  numerical  ex peri mentation  that  computed 
variation  of  tt  with  0m  (Figure  2)  seems  to  match  experimental  re¬ 
sults  most  faithfully  when  we  use  a  *  o*  *  1/2.  Effects  of  depar¬ 
tures  from  this  pair  of  values  are  so  pronounced,  in  fact,  that  our 
computations  seem  to  indicate  a  *  a*  =  1/2  represents  a  saddle  point 
in  closure-coefficient  space!  Thus,  we  conclude  that  a  and  a*  are 
equal  and  the  most  appropriate  value  is  1/2. 
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4.  SUBLAYER  ANALYSIS 


In  order  to  -facilitate  integration  of  the  eodel  equations  through 
the  viscous  sublayer ,  Me  oust,  at  a  ainiaue,  add  eolecular  diffu¬ 
sion  teres  to  the  equations  of  notion.  Potentially,  we  eight  al¬ 
so  have  to  alloe  the  various  closure  coefficients  to  be  functions 
of  viscosity  as  Mali.  In  this  section,  mm  use  perturbation  Meth¬ 
ods  to  analyze  viscous  sublayer  structure  predicted  by  the  new 
eodel,  including  effects  of  surface  roughness  and  surface  ease 
injection.  Note  that  mm  confine  our  analysis  to  the  mm  nodel  as 
results  of  Section  3  indicate  it  is  superior  to  the  other  nodels 
considered. 


4.1  PERTURBATION  SOLUTION 

Considering  the  constant-pressure  case  only,  convective  teres  are 
negligible  in  the  sublayer |  thus,  the  equations  of  notion  for  the 
new  nodel  (with  eolecular  diffusion  included)  sieplify  to* 


<  v+Vjj,  )  du/dy 

-  u2 

T 

(24) 

vT<du/dy> 

2  —  8*k  w  ♦ 

(d/dy)  [< v+o«vT>dk/dy]  ■ 

0 

(29) 

y (du/dy ) 

2  -  8  u>2  ♦ 

(d/dy> [ ( v+o  v^,)  du/dy]  « 

O 

(26) 

Five  boundary 

conditions  i 

re  needed,  two  of  which 

follow 

froe 

Matching  to  the  law  of  the  wall  as  y *  ->  •,  viz, 

Ic  ->  u2//|T*  and  u»  ->  uT//5*Ky  as  y*  ->  •  (27) 

Two  eore  boundary  conditions  follow  froe  "no  slip"  at  the  surface 
Mhich  inplies  that  u  and  k  vanish  at  y  »  0.  Thus, 


*  k  ■  0  at 


(2B) 


The  final  condition  is  sinil ar  to  that  deduced  in  earlier  studies 
<Refs  7-4)  where  we  have  found  that  for  perfect 1 y-eeooth  surf aces 
eolecular  diffusion  and  dissipation  balance  in  Equation  (2d)  and 
this  leads  to* 


in  ->  dv/8y2 


♦ 

y  ->  0 


(24) 


Morm  general  boundary  conditions  for  rough  surfaces  and  for  sur¬ 
faces  with  mass  injection  will  be  devised  in  Subsections  4. 2-4. 3. 
For  now,  we  focus  on  the  perf ectl y— smooth  surface. 

As  part  of  the  solution  to  Equations  (24-29),  we  obtain  the  con¬ 
stant  in  the  law  of  the  wall,  B,  where 

u+  ->  -  i,ny+  +  B  as  y+  ->  ®  (30) 

K 

We  solve  the  sublayer  equations  by:  (a)  adding  unsteady  terms  to 
the  right-hand  sides  of  Equations  (23-26);  (b)  making  an  initial 

guess  of  the  solution;  and  (c)  using  an  implicit,  time— marching, 
second-order — accurate  program  to  generate  the  long-time  solution 
in  which  the  unsteady  terms  tend  to  zero.  The  velocity  is  com¬ 
puted  at  each  timestep  using  the  fourth-order  Runge— Kutta  method. 
The  program  used  is  an  improved  version  of  that  developed  in  the 
study  by  Wilcox  and  Traci  (Ref  8). 

Using  this  program,  we  find  that  Equations  (24-29)  predict  that 
the  smooth-wall  value  of  B  is 


The  fact  that  this  value  is  well  within  the  scatter  of  measured 
values  of  B  strongly  suggests  that  no  further  viscous  modifica¬ 
tions  are  needed  for  this  model. 

Figures  5  and  6  compare  computed  and  measured  (Refs  2,16,17)  sub¬ 
layer  velocity  profiles  in  linear  and  semi log  coordinates.  The 
two  figures  show  that  computed  velocities  generally  fall  within 
experimental  data  scatter.  In  Figure  7,  we  compare  computed  and 
measured  turbulence  production  and  dissipation  terms.  Again,  pre¬ 
dictions  fall  well  within  experimental  error  bounds. 

Perhaps  the  only  deficiency  of  the  predicted  smooth— surf ace  sub¬ 
layer  structure  is  that  very  near  the  surface  the  model  predicts 

k  -v  y3*23  as  y  ->  0  (32) 


By  contrast,  the  Wi lcox— Rubesin  model  predicts  that  k  %  y1*  which 
suggests  that  k  %  <v,2>,  a  point  this  researcher  has  made  as  a 
more  physically  plausible  interpretation  than  saying  k  is  the  ki¬ 
netic  energy  of  the  turbulence.  By  letting  the  closure  coeffi¬ 
cient  8*  increase  as  a  function  of  turbulent  Reynolds  number,  Re^> 
defined  as  k/u>v,  it  is  possible  to  force  k  yH  but  then  we  find 
that  to  recover  B  ■  3,  at  least  two  other  closure  coefficients 
must  vary  with  Reip  (see  Appendix  B) .  Such  additional  complex  ity 
is  pointless  in  light  of  Figures  3  through  7. 
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Comparison  of  computed  and  measured  near-surface 
velocity  profiles  for  a  boundary  layer  over  a 
perfectly-smooth  surface. 


4.2  ROUGH-WALL  ANALYSIS 


A  key  advantage  of  the  k-oi2  and  k-cu  formulations  over  the  k— e 
■Formulation  is  the  -fact  that  u-oriented  equations  possess  solu¬ 
tions  in  which  the  value  of  u  may  be  arbitrarily  specified  at  the 
surface  (Ref  IS).  This  is  an  advantage  because  it  provides  a 
natural  way.  to  incorporate  effects  of  surface  roughness  through 
surface  boundary  conditions.  This  feature  of  the  equations  was 
originally  recognized  by  Saffman  (Ref  18). 


If  we  now  rewrite  the  surface  boundary  condition  (Equation  29)  on 
a)  as  follows. 


U) 


at  y  =  0 


(33) 


we  can  generate  sublayer  solutions  for  arbitrary  Sp  ,  including 

the  limiting  cases  S  ->  0  and  S  ->  ».  R 

A  R 

Figure  8  shows  the  computed  value  of  B  for  a  wide  range  of  values 

of  S  .  As  shown,  in  the  limit  Sr  ->  °°,  B  tends  to  S.  1.  To  de¬ 

termine  the  limiting  value  of  B  as  Sr  becomes  small,  we  replot 
the  numerical  results  on  semi log  paper  in  Figure  9.  As  shown,  an 
excellent  correlation  of  the  numerical  predictions  is  given  byt 

B  ->  8.4  +  -  jln^/lOO)  as  S0  ->  O  (34) 

K  A  A 

By  experimental  means,  Nikuradse  (Ref  19)  has  found  that  for  flow 
over  very  rough  surfaces; 

B  ->  8.3  +  i  Jln(l/kR)  ;  kj  =  uxkR/v  <35) 


where  kR  is  the  average  height  of  sand-grain  roughness  elements. 
(Note  that  in  our  computations  we  use  <  ■  0.41  while  Nikuradse 
used  k  *  0.40.)  Thus,  if  we  make  the  correlation! 


SR  -  100/kR  |  kR  »  1  (36) 


then  Equations  (34)  and  (33)  are  nearly  identical.  Figure  10  com¬ 
pares  computed  velocity  profiles  with  the  analytical  fit  obtained 
by  using  Equations  (34-33)  in  the  law  of  the  wall,  viz, 

u +  j  in(y/kp)  ♦  8.4  (37) 
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Figure  9*  Asymptotic  behavior  of  the  computed  constant  in 
the  law  of  the  wall,  B,  for  small  SB. 


Figure  10.  Computed  velocity  profiles  for  "completely-rough"  surfaces 


-for  three  of  the  computations.  The  correlation  is  nearly  exact. 
The  most  remarkable  fact  about  this  correlation  is  that  Equation 
(37)  is  the  form  the  law  of  the  wall  assumes  for  flow  over  "com- 
pletely-rough"  surfaces,  including  the  value  of  the  additive  con¬ 
stant  <8.4  and  8.5  differ  by  one  percent). 

By  making  a  qualitative  argument  based  on  flow  over  a  wavy  wall, 
Wilcox  and  Chambers  (Ref  20)  have  shown  that  for  small  roughness 
heights,  we  should  expect  to  have 


SR  *  <1/^,’ 


k+  ->  o 


(38) 


Comparison  with  Nikuradse’s  data  (Figure  11)  shows  that  the  fol¬ 
lowing  correlation  between  SR  and  kR  will  reproduce  measured 

effects  of  sand-grain  roughness  for  values  of  k+  up  to  about  400. 

n 


!(50/kR)2 

100/kJ 


kR  <  25 

kR  >  2S 


(39) 


4.3  EFFECTS  OF  SURFACE  MASS  INJECTION 

For  boundary  layers  with  surface  mass  injection,  the  introduction 
of  an  additional  velocity  scale  (v  *  normal  flow  velocity  at  the 
surface)  suggests  that  the  scalingWfor  u  at  the  surface  may  dif¬ 
fer  from  Equation  (33).  Andersen,  et  al  (Ref  17)  provide  further 
evidence  that  the  dissipation— r ate  boundary  condition  must  be  re¬ 
vised  when  mass  injection  is  present  by  showing,  from  correlation 
of  their  experimental  data,  that  both  <  and  B  are  functions  of 
v+  =*  vw/u  .  Because  our  rough-surf  ace  computations  of  the  pre¬ 
cedi  ngw subsection  show  that  the  value  of  B  is  strongly  affected 
by  the  surface  value  of  the  dissipation  rate,  this  suggests  that 
the  surface  value  of  u  will  depend  in  some  manner  upon  v+  .  Ex¬ 
amination  of  the  limiting  form  of  the  model  equations  for  y+->  00 
(i.e. ,  the  “wall-layer". . .Subsection  2.2)  shows  immediately  that 
the  effective  Karman  “constant”  varies  with  v*  according  tot 

W 

K  ■  </  <  1  5  v*)  (40) 


where  5  is  given  by 


■  3.11  +  0.61  £,ny+ 


(41) 
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1 


The  variation  of  ic  predicted  in  Equations  (40—41)  is  consistent 
with  the  Andersen,  et  al  data.  Including  appropriate  convective 
terns  in  Equations  (24-26),  we  have  performed  sublayer  computa¬ 
tions  for  the  cases  experimentally  documented  by  Andersen,  et  al . 
In  each  case,  the  surface  value  of  u)  has  been  given  by 

a)  =  - -  S  _  at  y  *  0  (42) 

V  D 

and  the  value  of  Sg  has  been  varied  until  optimum  agreement  be¬ 
tween  measured  and  computed  velocity  profiles  is  achieved.  The 
final  correlation  between  Sg  and  v+  is  shown  in  Figure  12.  The 
correlation  is  given  in  analytical  form  asi 


S 


B 


20 

VW(1+5VW> 


(43) 


Figure  13  displays  the  level  of  agreement  between  theory  and  ex¬ 
periment  using  Equations  (42)  and  (43). 


This  concludes  our  formulation  of  the  new  turbulence  model  and 
attending  surface  boundary  conditions.  In  the  following  section, 
we  focus  on  application  of  the  new  model  in  several  numerical 
computations. 
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S.  BOUNDARY-LAYER  APPLICATIONS 


We  turn  now  to  application  of  the  new  model  equations  to  a  total  of 
five  incompressible  boundary  layers.  First,  we  present  details  of 
the  manner  in  which  we  determine  boundary  conditions  for  the  turbu¬ 
lence-model  parameters  appropriate  at  the  boundary— 1 ayer  edge.  Then 
we  apply  the  model  to  the  constant-pressure  (flat-plate)  case.  The 
next  two  applications  are  to  boundary  layers  in  an  adverse  pressure 
gradient.  The  final  two  applications  are  for  boundary  layers  with 
surface  mass  transfer.  All  computations  have  been  done  with  the 
model  as  formulated  in  Sections  2  through  4  and,  with  the  viscous 
modifications  of  Appendix  B  included.  Graphical  results  are  almost 
identical  for  each  corresponding  pair  of  computations.  Thus,  only 
one  curve  is  presented  for  the  flow  properties  displayed. 


S. 1  BOUNDARY-LAYER-EDGE  CONDITIONS 

While  the  preceding  section  discussed  surface  boundary  conditions, 
we  have  not  yet  commented  on  boundary  conditions  appropriate  in  the 
freestream,  or  in  the  case  of  a  boundary  layer,  at  the  edge  of  the 
layer.  In  prior  computations  with  our  two-dimensional  boundary- 
layer  program,  EDDYBL  (Ref  21),  we  have  held  the  ratio  of  turbulent 
energy  to  the  square  of  the  mean  edge  velocity  constant  at  the  edge 
of  the  layer.  Additionally,  we  have  held  the  ratio  of  the  turbulent 
length  scale,  &  *  /k/u)  ,  to  the  local  boundary-1  ayer  thickness  con¬ 
stant.  These  two  boundary  conditions  are  more -or — less  consistent 
with  our  defect-layer  analysis.  However,  in  past  studies,  we  have 
found  that  these  boundary  conditions  can  lead  to  numerical  diffi¬ 
culties,  particularly  at  high  Reynolds  numbers.  In  order  to  avoid 
these  numerical  difficulties  in  this  study,  we  have  adopted  a  dif¬ 
ferent  approach. 

Ideally,  we  would  like  to  implement  "zero-gradient "  boundary  condi¬ 
tions  at  the  boundary-1 ayer  edge.  While  such  conditions  are  "clean" 
from  a  theoretical  point  of  view,  they  are  undesirable  from  a  nu¬ 
merical  point  of  view.  That  is,  the  conditions  we  have  used  in  past 
applications  are  of  the  Neumann  type  while  "zero-gradient "  condi¬ 
tions  are  of  the  Dirichlet  type.  Almost  universally,  convergence  of 
iterative  numerical  schemes  (EDDYBL  uses  an  iterative  scheme)  is 
much  slower  with  Dirichlet  conditions  than  with  Neumann  conditions. 

In  order  to  resolve  this  apparent  dilema,  we  appeal  directly  to  the 
equations  of  motion.  Beyond  the  boundar y-1 ayer  edge,  we  have  van¬ 
ishing  normal  gradients  so  that  the  equations  for  k  and  u  simplify 
to  the  following. 
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where  x  is  freestream  -flow  direction,  subscript  e  denotes  the  value 
at  the  boundary-layer  edge,  and  g  and  g*  are  defined  in  Equations 
(8).  The  solution  to  Equations  (44)  can  be  obtained  by  siaple  quad¬ 
rature,  independent  of  integrating  the  equations  of  eotion  through 
the  boundary  layer.  Once  ke  and  oi  are  determined  from  Equations 
(44),  it  is  then  possible  to  specify  Neumann-type  boundary  condi¬ 
tions  which  guarantee  zero  normal  gradients. 

Using  boundary  conditions  based  on  quadrature  of  Equations  (44),  we 
have  performed  a  series  of  numerical  experiments  to  test  solution 
sensitivity  to  the  boundary  conditions.  In  all  of  our  tests,  very 
little  solution  sensitivity  is  evident.  Most  importantly,  the  two 
turbulence  parameters  k  and  w  tend  smoothly  to  their  freestream 
values.  By  contrast,  the  boundary  conditions  we  have  used  in  past 
studies  produce  sharp  turbulent-nonturbulent  interfaces  which  some¬ 
times  hamper  solution  convergence.  The  new  boundary  conditions  thus 
appear  to  be  more  satisf actory. 


5.2  FLAT-PLATE  BOUNDARY  LAYER 

Our  first  application  is  the  const ant -pressure  (flat -pi ate)  boun¬ 
dary  layer.  While  this  application  does  not  provide  a  severe  test 
of  the  new  model,  it  is  nevertheless  necessary  to  be  sure  the  boun¬ 
dary-layer  program  has  been  coded  properly.  Additionally,  the  new 
model  wouldn't  be  of  much  use  as  a  predictive  tool  if'  it  were  inac¬ 
curate  for  this  simplest  of  all  boundary  layers. 

i 

The  computation  begins  at  a  plate-length  Reynolds  number,  Rex  ,  of 
one  million  and  continues  to  an  Rex  of  10.9  million.  A  total  of 
244  steps  are  taken  in  the  streamwise  direction.  Grid  points  normal 
to  the  surface  are  spaced  in  a  geometric  progression  with  a  grading 
ratio  of  111)  the  total  number  of  points  used  to  resolve  the  layer 
increases  from  61  initially  to  66  by  the  end  of  the  computation  (to 
accomodate  boundary-1 ayer  growth).  Numerical  experimentation  with 
the  current  version  of  EDDYBL  (Ref  22)  has  shown  that  our  solutions 
are  grid  independent  for  60  mesh  points.  Accuracy  to  within  3  per — 
cent  can  be  achieved  with  as  few  as  40  mesh  points. 

Figures  14  and  15  compare  computed  and  measured  skin  friction  and 
velocity  profiles,  respectively.  As  shown  in  Figure  14,  computed 
skin  friction  virtually  duplicates  correspond! ng  measurements  for 
the  entire  range  of  Reynolds  numbers  considered.  Figure  15  shows 
that  differences  between  computed  and  measured  velocity  profiles 
are  no  more  than  3  percent  of  scale  for  the  three  Reynolds  numbers 
indicated. 

Thus,  as  no  great  surprise,  the  new  model  is  quite  accurate  for  the 
flat-plate  boundary  layer.  In  the  next  two  applications,  we  address 
boundary  layers  experiencing  an  adverse  pressure  gradient. 
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Figure  15.  Comparison  of  computed  and  measured  velocity  profiles  for  a  flat-plate 
boundary  layer. 


S.3  BOUNDARY  LAYERS  WITH  ADVERSE  PRESSURE  GRADIENT 


We  consider  two  boundary  layers  with  adverse  pressure  gradient. 
The  first  case  has  a  moderate  adverse  pressure  gradient ,  the  ex¬ 
perimental  data  being  those  of  Bradshaw  (Ref  2... Case  3300).  The 
second  case  has  increasingly  adverse  pressure  gradient,  the  exper — 
i mental  data  being  those  of  Samuel  and  Joubert  (Ref  10. ..Flow  S2) . 

For  the  Bradshaw  case,  a  total  of  90  streamwise  steps  are  needed 
to  march  from  x  *  2.5  ft  to  x  *  7.0  ft,  corresponding  to  Rex  in¬ 
creasing  from  about  2  million  to  about  4  million.  A  geometric- 
progression  ratio  of  1SH  is  used  to  construct  the  grid  normal  to 
the  surface  with  a  total  of  61  points  throughout  the  computation. 

Figures  16  and  17  compare  computed  and  measured  skin  friction  and 
a  velocity  profile.  Inspection  of  both  graphs  shows  that  differ — 
ences  between  theory  and  experiment  nowhere  exceed  5  percent  for 
this  flow. 

In  the  Samuel -Joubert  case,  we  use  100  streamwise  steps  to  march 
from  x  —  1  a  to  xs  3.40  m,  corresponding  to  an  Rex  range  of  about 
2  million  to  4  million.  To  achieve  an  extremely  accurate  numer — 
ical  solution  in  this  case,  we  have  used  a  geometric— progression 
ratio  of  7. 87  and  a  total  of  81  mesh  points.  Results  differ  from 
a  61-point  computation  by  less  than  a  half  a  percent. 

Figures  18  and  19  compare  computed  and  measured  skin  friction  and 
two  velocity  profiles  for  this  flow.  Computed  and  measured  skin 
friction  differ  by  less  than  7  percent  of  scale.  The  velocity  pro¬ 
files  at  x  »  2.87  m  are  within  5  percent  while  those  at  x  -  3.40  m 
differ  by  no  more  than  9  percent.  Although  our  predictions  for 
this  case  differ  from  measurements  a  bit  more  than  the  other  two 
cases,  note  that  our  predictions  for  this  flow  at  the  Second  Stan¬ 
ford  Olympiad  typically  differed  from  measurements  by  as  much  as 
25  percent,  and  our  prediction  was  one  of  the  best  submitted. 


5.4  BOUNDARY  LAYERS  WITH  SURFACE  NASS  TRANSFER 

As  our  final  two  applications  of  the  model,  we  consider  two  cases 
with  surface  mass  transfer.  Both  cases  were  included  in  the  1981 
Stanford  Olympics  (Ref  10)  and  data  for  both  cases  were  taken  by 
Andersen,  et  al  (Ref  17).  The  first  case  has  surface  mass  injec¬ 
tion  rate,  v  .  given  by  .00375  U@.  The  second  case. has  surface 
mass  removal  Vate,  v ^  given  by  -.00375  Ue.  For  the  blowing  case 
pressure  is  constant,  while  the  suction  case  has  a  slight  adverse 
pressure  gradient. 

Figures  20  and  21  compare  computed  and  measured  skin  friction  and 
velocity  profiles,  respectively,  for  the  blowing  case.  As  shown, 
computed  and  measured  skin  friction  differ  by  less  than  4  percent 
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Figure  18.  Comparison  of  computed  and  measured  skin  friction 
for  the  Samuel-Joubert  adverse-pressure-gradient 
boundary  layer. 
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Figure  20.  Comparison  of  computed  and  measured  skin  friction 
for  a  boundary  layer  with  constant  pressure  and 
uniform  mass  injection;  v  ■  .00375  U  . 


Figure  21.  Comparison  of  computed  and  measured  velocity  profiles  for  a  boundary  layer 
with  constant  pressure  and  uniform  mass  Injection;  v  =  .00375  U  . 


of  meals  while  computed  and  measured  velocity  profiles  are  within 
3  pmreant  of  mach  othmr. 

Figurn  22  and  23  compare  computed  and  measured  skin  friction  and 
velocity  profiles,  raspocti vsly,  for  the  suction  cats.  As  shown, 
aftsr  tho  initial  transient  (caused  by  a  poor  choice  of  initial 
k  and  profiles),  computed  skin  friction  asymptotes  to  a  value  4 
percent  higher  than  measured.  The  computed  velocity  profile  lies 
within  3  percent  of  corresponding  measurements. 


l 
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Figure  22.  Comparison  of  computed  and  measured  skin  friction 
for  a  boundary  layer  with  mild  adverse  pressure 
gradient  and  uniform  suction;  v  ■  -.00375  U  . 


Figure  23.  Comparison  of  computed  and  measured  velocity  profiles  for  a  boundary  layer  with 
mild  adverse  pressure  gradient  and  uniform  suction;  v  =  -.00375  U  . 


6.  SUMMARY  AND  CONCLUSIONS 


The  primary  objectives  of  this  study  have  been  accomplished,  viz, 
we  have  made  a  critical  review  of  closure  approximations  used  in 
two— equation  turbulence  models  and  determined  what  appears  to  be 
an  optimum  choice  of  dependent  variables.  As  a  result,  we  have 
developed  a  new  two-equation  model  which  shows  promise  of  being 
more  accurate  for  boundary  layers  in  an  adverse  pressure  gradient 
than  any  other  similar  model. 

As  in  our  prior  turbulence  modeling  efforts,  we  have  made  exten¬ 
sive  use  of  perturbation  methods  (Sections  3  and  4).  In  contrast 
to  prior  studies,  our  analysis  of  the  defect  layer  includes  pres¬ 
sure  gradient.  As  discussed  in  Section  3,  limiting  the  defect- 
layer  analysis  to  the  constant-pressure  case  displays  little  dif¬ 
ference  amongst  the  various  two-equation  models  in  general  usage. 
However,  as  soon  as  an  adverse  pressure  gradient  is  included,  the 
models  exhibit  large  differences.  As  a  general  observation,  the 
Second  Stanford  Olympiad  demonstrated  that  modern  turbulence  mo¬ 
dels  are  not  much  more  accurate  than  those  in  use  13  years  ago  if 
the  flow  of  concern  is  a  boundary  layer  in  adverse  pressure  gra¬ 
dient.  The  analysis  of  Section  3  indicates  why  this  is  true  and, 
with  the  introduction  of  the  new  model,  offers  the  basis  for  de¬ 
velopment  of  new  models  which  are  accurate  for  such  flows. 

The  model  thus  far  has  been  tested  for  five  incompressible  boun¬ 
dary  layers.  Certainly  much  more  testing  is  needed  before  the 
model  can  be  offered  for  general  use.  Additional  testing  will  be 
done  in  future  research  efforts.  However,  additional  development 
of  the  model  will  be  needed  before  such  tests  can  or  should  be 
done.  Specifically,  the  model  must  be  generalized  for  compres¬ 
sibility  and,  even  more  importantly,  the  constitutive  relation 
between  Reynolds  stress  and  mean-flow  properties  must  be  revised. 
The  relation  proposed  in  Equation  (6)  fails  to  predict  anisotropy 
of  the  normal  stresses,  and  does  not  account  for  streamline  cur¬ 
vature  effects.  Additionally,  it  is  not  at  all  clear  that  model 
predictions  will  bpar  any  relation  to  physical  reality  for  flows 
which  are  unsteady.  Because  unsteady  boundary  layers  are  of  key 
interest  in  the  overall  Contract  objectives,  future  development 
efforts  must  include  unsteady  flow  analysis  and  application. 


APPENDIX  A:  DEFECT-LAYER  EQUATIONS 


In  this  Appendix,  we  present  details  of  the  formal  perturbation  ex¬ 
pansion  solution  to  the  defect-layer  equations.  First,  we  summarize 
the  three  turbulence  models  under  consideration.  Next,  we  outline 
the  form  of  the  perturbation  expansions  and  state  the  equations  for 
the  leading  order  terms  in  the  expansions.  Then,  we  present  boun¬ 
dary  conditions  used  in  solving  the  defect— layer  equations.  Final¬ 
ly,  we  indicate  the  manner  in  which  skin  friction  and  wake  strength 
can  be  extracted  from  the  defect-layer  solution. 


A. 1  TURBULENCE  MODELS  UNDER  CONSIDERATION 

In  analyzing  the  defect  layer,  we  focus  on  three  turbulence  models, 
viz,  the  new  model  postulated  in  Section  2,  the  Ni lcox-Rubesin  mo¬ 
del  (Ref  9)  and  the  Jones— Launder  model  (Ref  5).  For  all  three  of 
the  models  we  must  solve  the  equations  for  mean  mass  and  momentum 
conservation,  an  equation  for  a  turbulent  energy  scale  and  an  equa¬ 
tion  for  a  turbulent  dissipation  scale.  For  all  three  models,  the 
first  three  equations  assume  the  following  form. 
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Note  that  Equations  (A. 1-A.3)  do  not  include  molecular  viscosity. 
This  is  a  valid  approximation  in  the  defect  layer  as  the  eddy  vis¬ 
cosity  is  proportional  to  Ue6*,  where  U@  is  the  boundary-1 ayer — edge 
velocity  and  5*  is  displacement  thickness.  Hence,  the  ratio  of  mo¬ 
lecular  to  eddy  viscosity  varies  inversely  with  displacement-thick¬ 
ness  Reynolds  number  and  is  thus  very  small.  The  difference  amongst 
the  three  models  is  the  way  the  dissipation,  e  ,  and  the  eddy  vis¬ 
cosity,  v,p,  are  computed.  , 


A.  1 . 1  New  Model 

For  this  model,  in  addition  to  Equations  (A. 1-A.3)  we  havei 
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Leaving  o  and  a*  undetermined,  the  values  for  the  remaining  closure 
coefficients  arei 
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A. 1.2  Wilcox-Rubesin  Model 

In  this  model,  the  additional  equations  are  as  foil 
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A. 1 . 3  Jones -Launder  Model 

In  this  model  tm  compute  dissipation,  e  , 
additional  equations  arei 
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Values  for  the  closure  coefficients  arei 
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A. 2  EXPANSION  PROCEDURE 


Following  the  formulation  of  Milcox  and  Traci  (Ref  8),  we  introduce 
a  streaef unction  i|>  ,  and  seek  a  perturbation  solution  of  the  form 
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Inserting  Equations  (A.8-A.12)  into  Equations  (A. 1— A. 7) ,  neglecting 
higher — order  teres,  letting  N  denote  diaensionless  eddy  viscosity, 
and  defining 
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we  obtain  the  following  equations. 
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The  final  equation  for  each  eodel  is  differenti  the  equations  arei 
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Ni lcox-Rubesin  Model 
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Jones-Launder  Model 
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where  the  parameters  ctrp.  Sip,  Oip  and  u>p  are  defined  in  terms  of  6*, 
uT  and  skin  friction,  Cpj  *  2u*/u|,  i.e. , 
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Equations  (A.  14— A.  18)  will  have  self-similar  solutions  only  if  a<p, 
S>p  and  uiip  are  independent  of  x.  As  noted  by  Bush  and  Fendell 
(Ref  14),  for  Re^v>>  1,  uT  varies  sufficiently  slowly  that  we  have 
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and,  in  addition,  the  shape  factor  to  leading  order  approaches  one 
which  (from  inspection  of  the  momentum-integral  equation)  implies 
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Thus,  self-similar  solutions  exist  provided  the  paramet 
dependent  of  x. 
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In  summary,  the  defect-layer  equations  for  the  leading  order  ten 
in  the  perturbation  expansions  becomes 
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A. 3  BOUNDARY  CONDITIONS 

At  the  outer  edge  of  the  defect  layer,  we  pose  the  conditions  that 
the  velocity  equal  the  freestream  velocity.  Additionally,  we  let 
the  turbulent  energy  assume  a  small  value  and  insist  that  the  tur¬ 
bulent  length  scale  have  zero  slope  at  the  boundary-1 ayer  edge.  In 
their  defect— layer  analysis,  Wilcox  and  Traci  used  these  boundary 
conditions  as  well  as  explicitly  prescribing  both  1^  and  WQ.  Thus, 

U1  *  0  ) 

K0  ■  small  value  >  at  n  •  ng  (A. 25) 
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Approaching  the  surface,  we  eust  formally  match  to  the  law  of  the 
wall.  Matching  is  a  bit  different  for  each  model  but  is  neverthe¬ 
less  straightf orwardj  details  of  the  algebra  will  thus  be  omitted 
in  the  interest  of  brevity.  The  limiting  forms  we  have  used  for 
n  ->  0  follow. 
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Additionally,  the  coefficient  u  is  determined  from  the  integral 
constraint  for  mass  conservation?  viz. 


(A.  30) 


To  implement  Equation  (A. 30)  in  the  numerical  solution  of  the  de¬ 
fect-layer  equations,  we  proceed  as  follows.  Because  the  velocity 
is  singular  at  n  ■  0,  we  integrate  the  numerical  velocity  profile 
from  a  point  above  the  surface  which  we  denote  by  (typically  of 
order  . OOl  to  .01)  to  the  edge  of  the  boundary  layer,  He •  Then, 
we  integrate  the  asymptotic  profile  for  lij^  (ri)  given  in  Equations 
(A.  26)  from  n  *  0  to  n  ■  r^.  The 'latter  integration  involves 
the  unknown  coefficient  uQ .  Finally,  the  sum  of  these  two  inte¬ 
grals  must  be  unity  by  virtue  of  Equation  (A. 30). 


A. 4  SKIN  FRICTION  AND  WAKE  STRENGTH 

It  is  possible  to  determine  the  skin  friction  implied  by  the  solu¬ 
tion  to  the  defect— layer  equations  by  formal  matching  to  the  sub¬ 
layer  asymptotic  velocity  profile.  Considering  only  leading  order 
terms,  this  means  we  sayt 
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Finally,  combining  Coles'  composite  profile 
Equation  (A. 32)  and  evaluating  the  resulting 
d ary-layer  edge,  we  arrive  at  the  following 
strength,  } i 


[Equation  (19)3  with 
equation  at  the  boun- 
ex press ion  for  wake 


*  ■  Huo  -  ln”J 


(A.  33) 
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APPENDIX  Bs  VISCOUS  MODIFICATIONS 


As  noted  in  Section  4,  very  close  to  the  surface  the  nett  model 
equations  predict  the  turbulent  energy  tends  to  zero  as  distance 
to  the  3.23  power.  By  contrast,  if  the  quantity  k  is  turbulent 
kinetic  energy  we  would  expect  to  have  k  ^  y2  as  y  ->  0.  Al¬ 
ternatively,  if  k  is  kinetic  energy  of  surf ace-noreal  fluctua¬ 
tions,  we  would  expect  to  have  k  ^  y4  as  y  ->  O.  Either  way, 
we  do  not  expect  to  have  a  non-integer  exponent  such  as  3.23. 

While,  on  the  one  hand,  this  feature  of  model  predictions  is  of 
minor  interest,  there  are  other  shortcomings  of  the  model  as  for¬ 
mulated  in  Sections  2—4  which  are  not  necessarily  minor.  For  ex¬ 
ample,  if  transition  from  laminar  to  turbulent  flow  is  important 
for  a  given  application,  the  model  must  be  altered  in  order  to 
achieve  physically  realistic  predictions.  Addi tionally,  predic¬ 
tions  for  1 ow-Reynol ds-number  turbulent  flows  conceivably  are 
inaccurate  if  no  provision  is  made  for  viscous  effects  upon  the 
various  closure  coefficients.  To  make  provision  for  extending  the 
model’s  applicability  to  such  flows,  this  Appendix  summarizes  the 
required  viscous  modifications  to  the  model  equations. 
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and  the  closure  coefficients  aret 


8  -  3/40,  a  -  1/2,  ct*  -  1/2 

0*  "  100  t1  +  f  exp(-ReT/4)] 
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Additionally,  the  quantities  Re^ 


* 


X  ,  R,  and  R  are  given  by 

k  m 


Ratp  -  k/uv  j  X  »  /6/59  ;  RJc  *  1  ;  R^  *  9/4  (B.S) 


The  rationale  -for  the  modified  closure  coefficients  follows. 

Considering  first  the  coefficient  6*  ,  the  model -predi cted  very 
near  wall  behavior  of  k  now  becomes  k  %  y4  as  y  — >  0.  Thus,  our 
prediction  indicates  k  is  most  appropriately  identified  as  being 
the  kinetic  energy  of  the  fluctuations  normal  to  the  surface.  In 
addition,  the  coefficient  4  in  the  exponential  damping  term  has 
been  chosen  to  yield  satisfactory  agreement  with  the  experimental 
data  of  Batchelor  and  Townsend  (Ref  23)  for  the  late— ti me  decay 
of  homogeneous  turbulence. 

Turning  now  to  the  coefficients  Y  and  Y*  ,  Wilcox  and  Rubesin 
(Ref  9)  found  similar  modifications  necessary  in  their  model  in 
order  to  achieve  physically  acceptable  predictions  in  the  viscous 
sublayer.  The  indicated  choice  of  the  coefficient  X  guarentees 
that  in  a  Blasius  boundary  layer,  turbulent  energy  will  only  be 
amplified  for  Reynolds  nuiber  exceeding  the  linear-stability  min¬ 
imum-critical  Reynolds  number.  The  values  for  R  and  R  have 
been  determined  by  analyzing  the  smooth-surface  viscous  sublayer. 
The  values  quoted  insure  the  model  predicts  the  additive  constant 
in  the  law  of  the  wall  is  9.0. 

We  have  rerun  the  Section  9  boundary-layer  test  cases  using  Equa¬ 
tions  (B. 1-B.9).  In  none  of  the  cases  have  we  detected  signifi¬ 
cant  differences  from  results  obtained  when  no  viscous  modifica¬ 
tions  appear  in  the  model  equations.  However,  we  do  notice  a  20% 
Increase  in  computing  time  because  of  the  many  exponentials  which 
must  be  evaluated  at  every  mesh  point.  We  thus  conclude  that  the 
viscous  modif i cat ions  proposed  in  Equations  (B. 1-B.9)  add  little 
to  the  model  for  ful ly-turbulent  boundary-1 ayer  applications. 
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